clip-02-05.jl

Load Julia packages (libraries) needed for the snippets in chapter 0

using StatisticalRethinking, Optim
gr(size=(600,300));

Plots.GRBackend()

snippet 3.2

p_grid = range(0, step=0.001, stop=1)
prior = ones(length(p_grid))
likelihood = [pdf(Binomial(9, p), 6) for p in p_grid]
posterior = likelihood .* prior
posterior = posterior / sum(posterior)
samples = sample(p_grid, Weights(posterior), length(p_grid));
samples[1:5]

5-element Array{Float64,1}:
 0.6  
 0.613
 0.388
 0.691
 0.548

snippet 3.3

Draw 10000 samples from this posterior distribution

N = 10000
samples = sample(p_grid, Weights(posterior), N);

10000-element Array{Float64,1}:
 0.565
 0.671
 0.591
 0.746
 0.736
 0.602
 0.728
 0.753
 0.741
 0.53 
 ⋮    
 0.403
 0.744
 0.373
 0.573
 0.626
 0.57 
 0.732
 0.47 
 0.662

In StatisticalRethinkingJulia samples will always be stored in an MCMCChains.Chains object.

chn = MCMCChains.Chains(reshape(samples, N, 1, 1), ["toss"]);

Object of type Chains, with data of type 10000×1×1 Array{Float64,3}

Iterations        = 1:10000
Thinning interval = 1
Chains            = 1
Samples per chain = 10000
parameters        = toss

parameters
      Mean    SD   Naive SE  MCSE     ESS  
toss 0.6366 0.1395   0.0014 0.0015 8606.187

Describe the chain

describe(chn)

Iterations        = 1:10000
Thinning interval = 1
Chains            = 1
Samples per chain = 10000
parameters        = toss

Empirical Posterior Estimates
───────────────────────────────────────────
parameters
      Mean    SD   Naive SE  MCSE     ESS
toss 0.6366 0.1395   0.0014 0.0015 8606.187

Quantiles
───────────────────────────────────────────
parameters
      2.5% 25.0% 50.0% 75.0% 97.5%
toss 0.347 0.543 0.646 0.741 0.877

Plot the chain

plot(chn)

snippet 3.4

Create a vector to hold the plots so we can later combine them

p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
p[1] = scatter(1:N, samples, markersize = 2, ylim=(0.0, 1.3), lab="Draws")